library(tidyverse)
library(readstata13)
library(foreign)
library(readxl)
library(writexl)
library(lubridate)
library(PanelMatch)
library(RColorBrewer)
library(stargazer)

# Set working directory
setwd("/Volumes/Backup Plus/Storage-OrganizedFiles/Article1-ReplicationMaterials_FullReplication/2.MainResults-Tables&Figures")
getwd()

# Charge the results of the main models
load("OMC-ResultsMainModel_Updated.RData")

################################################################
# Figure 1: Preparing the Data
################################################################

dat <- data_final %>%
  group_by(year) %>%
  summarise(Protest = sum(iv2.p_strat),
            Sanctuary = sum(iv3.s_strat),
            VSP = sum(iv4b.v_strat))

dat1 <- dat %>%
  gather(Strategy, Municipalities, Protest:VSP) %>%
  rename(Year = year)

my_colours = list(
  my_greys = c("#1A242F", "#474F58", "#757B82", "#A3A7AB", "#D1D3D5"),
  strategies = c("#377eb8", "#4daf4a","#e41a1c")
)

my_palettes = function(name, n, all_palettes = my_colours, type = c("discrete", "continuous")) {
  palette = all_palettes[[name]]
  if (missing(n)) {
    n = length(palette)
  }
  type = match.arg(type)
  out = switch(type,
               continuous = grDevices::colorRampPalette(palette)(n),
               discrete = palette[1:n]
  )
  structure(out, name = name, class = "palette")
}

my_palettes("strategies", type = "discrete")
my_palettes("my_greys", type = "continuous")

scale_colour_cvi_d = function(name) {
  ggplot2::scale_colour_manual(values = my_palettes(name,
                                                     type = "discrete"))
}

scale_fill_cvi_d = function(name) {
  ggplot2::scale_fill_manual(values = my_palettes(name,
                                                   type = "discrete"))
}

figure2 <- ggplot(data=dat1, aes(x=Year, y=Municipalities, fill=Strategy)) +
  geom_bar(stat = "identity", position=position_dodge(), width = 1) +
  scale_fill_cvi_d("strategies")

pdf(file = "Figure1-Final.pdf",
    width = 6, 
    height = 4)

figure2 + theme_test() + theme(text = element_text(size = 10)) + theme(legend.position="bottom")

dev.off()

################################################################
# Figure 2: Preparing the Data
################################################################

x4 <- summary(pe.results5.p2_CBPSm5)
x5c <- summary(pe.results5.s2_maha5b)
x1b <- summary(pe.results5.v2_CBPSw5)

data4a <- as.data.frame(x4$summary)
data4a <- data4a %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4a <- data4a %>%
  mutate(time = c(0, 1, 2, 3, 4))

data5c <- as.data.frame(x5c$summary)
data5c <- data5c %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5c <- data5c %>%
  mutate(time = c(0, 1, 2, 3, 4))

data1b <- as.data.frame(x1b$summary)
data1b <- data1b %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1b <- data1b %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Figure 2: Plot
################################################################

pdf(file = "Figure2-Final.pdf",
    width = 6, 
    height = 3)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4, 1.5, 0))
par(mfrow = c(1, 3))

plot(x = data4a$time,
     y = data4a$estimate, 
     main = "",
     ylim = c(-15, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-15, 15, 6),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data4a$time, x1 = data4a$time,
         y0 = data4a$lower_b, y1 = data4a$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5c$time,
     y = data5c$estimate, 
     main = "",
     ylim = c(-15, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-15, 15, 6),
     pch = 16, cex.axis=1.25)
segments(x0 = data5c$time, x1 = data5c$time,
         y0 = data5c$lower_b, y1 = data5c$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data1b$time,
     y = data1b$estimate, 
     main = "",
     ylim = c(-15, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-15, 15, 6),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1b$time, x1 = data1b$time,
         y0 = data1b$lower_b, y1 = data1b$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)


mtext("Protest", line = -1.8, at = 0.18, outer = TRUE, cex = 1)
mtext("Sanctuary", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("VSP", line = -1.8, at = 0.85, outer = TRUE, cex = 1)

mtext("Estimated Effect on Violence", side = 2, line = .3, at = 0.5, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years Relative to the Administration of Treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Figure 3: Preparing the Data
################################################################

x4d_1 <- summary(pe.results5.p2_CBPSm5_mod1[["Nondominant"]])
x5e_1 <- summary(pe.results5.s2_maha5b_mod1[["Nondominant"]])
x1f_1 <- summary(pe.results5.v2_CBPSw5_mod1[["Nondominant"]])
x4d_2 <- summary(pe.results5.p2_CBPSm5_mod1[["Dominant"]])
x5e_2 <- summary(pe.results5.s2_maha5b_mod1[["Dominant"]])
x1f_2 <- summary(pe.results5.v2_CBPSw5_mod1[["Dominant"]])

data4d1 <- as.data.frame(x4d_1$summary)
data4d1 <- data4d1 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4d1 <- data4d1 %>%
  mutate(time = c(0, 1, 2, 3, 4))

data5e1 <- as.data.frame(x5e_1$summary)
data5e1 <- data5e1 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5e1 <- data5e1 %>%
  mutate(time = c(0, 1, 2, 3, 4))

data1f1 <- as.data.frame(x1f_1$summary)
data1f1 <- data1f1 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1f1 <- data1f1 %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4d2 <- as.data.frame(x4d_2$summary)
data4d2 <- data4d2 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4d2 <- data4d2 %>%
  mutate(time = c(0, 1, 2, 3, 4))

data5e2 <- as.data.frame(x5e_2$summary)
data5e2 <- data5e2 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5e2 <- data5e2 %>%
  mutate(time = c(0, 1, 2, 3, 4))

data1f2 <- as.data.frame(x1f_2$summary)
data1f2 <- data1f2 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data1f2 <- data1f2 %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Figure 3: Plot
################################################################
pdf(file = "Figure3-Final.pdf",
    width = 6, 
    height = 6)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4, 1.5, 0))
par(mfrow = c(2, 3))

plot(x = data4d1$time,
     y = data4d1$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data4d1$time, x1 = data4d1$time,
         y0 = data4d1$lower_b, y1 = data4d1$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5e1$time,
     y = data5e1$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, cex.axis=1.25)
segments(x0 = data5e1$time, x1 = data5e1$time,
         y0 = data5e1$lower_b, y1 = data5e1$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data1f1$time,
     y = data1f1$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1f1$time, x1 = data1f1$time,
         y0 = data1f1$lower_b, y1 = data1f1$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)


plot(x = data4d2$time,
     y = data4d2$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data4d2$time, x1 = data4d2$time,
         y0 = data4d2$lower_b, y1 = data4d2$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5e2$time,
     y = data5e2$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     pch = 16, cex.axis=1.25)
segments(x0 = data5e2$time, x1 = data5e2$time,
         y0 = data5e2$lower_b, y1 = data5e2$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data1f2$time,
     y = data1f2$estimate, 
     main = "",
     ylim = c(-20, 20), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 20, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data1f2$time, x1 = data1f2$time,
         y0 = data1f2$lower_b, y1 = data1f2$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

mtext("Protest", line = -1.8, at = 0.18, outer = TRUE, cex = 1)
mtext("Sanctuary", line = -1.8, at = 0.51, outer = TRUE, cex = 1)
mtext("VSP", line = -1.8, at = 0.85, outer = TRUE, cex = 1)

mtext("Zone 4 = 0", side = 2, line = .3, at = 0.75, 
      outer = TRUE, cex = 1)
mtext("Zone 4 = 1", side = 2, line = .3, at = 0.25, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years Relative to the Administration of Treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Figure 4: Preparing the data
################################################################

x4g_1 <- summary(pe.results5.p2_CBPSm5_mod2[["Low"]])
x5h_1 <- summary(pe.results5.s2_maha5b_mod2[["Low"]])
x4g_2 <- summary(pe.results5.p2_CBPSm5_mod2[["High"]])
x5h_2 <- summary(pe.results5.s2_maha5b_mod2[["High"]])

data4g1 <- as.data.frame(x4g_1$summary)
data4g1 <- data4g1 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4g1 <- data4g1 %>%
  mutate(time = c(0, 1, 2, 3, 4))

data5h1 <- as.data.frame(x5h_1$summary)
data5h1 <- data5h1 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5h1 <- data5h1 %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4g2 <- as.data.frame(x4g_2$summary)
data4g2 <- data4g2 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4g2 <- data4g2 %>%
  mutate(time = c(0, 1, 2, 3, 4))

data5h2 <- as.data.frame(x5h_2$summary)
data5h2 <- data5h2 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5h2 <- data5h2 %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Figure 4: Plot
################################################################
pdf(file = "Figure4-Final.pdf",
    width = 5, 
    height = 6)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4, 1.5, 0))
par(mfrow = c(2, 2))

plot(x = data4g1$time,
     y = data4g1$estimate, 
     main = "",
     ylim = c(-25, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-25, 15, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data4g1$time, x1 = data4g1$time,
         y0 = data4g1$lower_b, y1 = data4g1$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5h1$time,
     y = data5h1$estimate, 
     main = "",
     ylim = c(-25, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-25, 15, 4),
     pch = 16, cex.axis=1.25)
segments(x0 = data5h1$time, x1 = data5h1$time,
         y0 = data5h1$lower_b, y1 = data5h1$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)


plot(x = data4g2$time,
     y = data4g2$estimate, 
     main = "",
     ylim = c(-25, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-25, 15, 4),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data4g2$time, x1 = data4g2$time,
         y0 = data4g2$lower_b, y1 = data4g2$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5h2$time,
     y = data5h2$estimate, 
     main = "",
     ylim = c(-25, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-25, 15, 4),
     pch = 16, cex.axis=1.25)
segments(x0 = data5h2$time, x1 = data5h2$time,
         y0 = data5h2$lower_b, y1 = data5h2$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)


mtext("Protest", line = -1.8, at = 0.26, outer = TRUE, cex = 1)
mtext("Sanctuary", line = -1.8, at = 0.77, outer = TRUE, cex = 1)

mtext("Vote for the Left: Low", side = 2, line = .3, at = 0.75, 
      outer = TRUE, cex = 1)
mtext("Vote for the Left: High", side = 2, line = .3, at = 0.25, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years Relative to the Administration of Treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Figure 5: Preparing the data
################################################################

x4i_1 <- summary(pe.results5.p2_CBPSm5_mod3[["No Recruitment"]])
x5j_1 <- summary(pe.results5.s2_maha5b_mod3[["No Recruitment"]])
x4i_2 <- summary(pe.results5.p2_CBPSm5_mod3[["Recruitment"]])
x5j_2 <- summary(pe.results5.s2_maha5b_mod3[["Recruitment"]])

data4i1 <- as.data.frame(x4i_1$summary)
data4i1 <- data4i1 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4i1 <- data4i1 %>%
  mutate(time = c(0, 1, 2, 3, 4))

data5j1 <- as.data.frame(x5j_1$summary)
data5j1 <- data5j1 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5j1 <- data5j1 %>%
  mutate(time = c(0, 1, 2, 3, 4))


data4i2 <- as.data.frame(x4i_2$summary)
data4i2 <- data4i2 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data4i2 <- data4i2 %>%
  mutate(time = c(0, 1, 2, 3, 4))

data5j2 <- as.data.frame(x5j_2$summary)
data5j2 <- data5j2 %>%
  rename(lower_b = "2.5%",
         upper_b = "97.5%")
data5j2 <- data5j2 %>%
  mutate(time = c(0, 1, 2, 3, 4))

################################################################
# Figure 5: Plot
################################################################

pdf(file = "Figure5-Final.pdf",
    width = 6, 
    height = 7)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4, 1.5, 0))
par(mfrow = c(2, 2))

plot(x = data4i1$time,
     y = data4i1$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 10, 6),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data4i1$time, x1 = data4i1$time,
         y0 = data4i1$lower_b, y1 = data4i1$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5j1$time,
     y = data5j1$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 10, 6),
     pch = 16, cex.axis=1.25)
segments(x0 = data5j1$time, x1 = data5j1$time,
         y0 = data5j1$lower_b, y1 = data5j1$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)


plot(x = data4i2$time,
     y = data4i2$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 10, 6),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data4i2$time, x1 = data4i2$time,
         y0 = data4i2$lower_b, y1 = data4i2$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5j2$time,
     y = data5j2$estimate, 
     main = "",
     ylim = c(-20, 10), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-20, 10, 6),
     pch = 16, cex.axis=1.25)
segments(x0 = data5j2$time, x1 = data5j2$time,
         y0 = data5j2$lower_b, y1 = data5j2$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)


mtext("Protest", line = -1.8, at = 0.26, outer = TRUE, cex = 1)
mtext("Sanctuary", line = -1.8, at = 0.77, outer = TRUE, cex = 1)

mtext("No Rebel Recruitment", side = 2, line = .3, at = 0.75, 
      outer = TRUE, cex = 1)
mtext("Rebel Recruitment", side = 2, line = .3, at = 0.25, 
      outer = TRUE, cex = 1)
mtext(1, text = "Years Relative to the Administration of Treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1)

dev.off()

################################################################
# Figure 4: Alternative Plot
################################################################

pdf(file = "Figure4-Alt.pdf",
    width = 8, 
    height = 7)

par(mar = c(1.5, 2, 2 , 1), oma = c(4, 4, 1.5, 0))
par(mfrow = c(2, 4))

plot(x = data4g1$time,
     y = data4g1$estimate, 
     main = "",
     ylim = c(-10, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-10, 15, 5),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data4g1$time, x1 = data4g1$time,
         y0 = data4g1$lower_b, y1 = data4g1$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4g2$time,
     y = data4g2$estimate, 
     main = "",
     ylim = c(-10, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-10, 15, 5),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data4g2$time, x1 = data4g2$time,
         y0 = data4g2$lower_b, y1 = data4g2$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4i1$time,
     y = data4i1$estimate, 
     main = "",
     ylim = c(-10, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-10, 15, 5),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data4i1$time, x1 = data4i1$time,
         y0 = data4i1$lower_b, y1 = data4i1$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data4i2$time,
     y = data4i2$estimate, 
     main = "",
     ylim = c(-10, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-10, 15, 5),
     xlab = "", ylab = "",
     pch = 16, cex.axis=1.25)
segments(x0 = data4i2$time, x1 = data4i2$time,
         y0 = data4i2$lower_b, y1 = data4i2$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)


plot(x = data5h1$time,
     y = data5h1$estimate, 
     main = "",
     ylim = c(-25, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-25, 15, 4),
     pch = 16, cex.axis=1.25)
segments(x0 = data5h1$time, x1 = data5h1$time,
         y0 = data5h1$lower_b, y1 = data5h1$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5h2$time,
     y = data5h2$estimate, 
     main = "",
     ylim = c(-25, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-25, 15, 4),
     pch = 16, cex.axis=1.25)
segments(x0 = data5h2$time, x1 = data5h2$time,
         y0 = data5h2$lower_b, y1 = data5h2$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5j1$time,
     y = data5j1$estimate, 
     main = "",
     ylim = c(-25, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-25, 15, 4),
     pch = 16, cex.axis=1.25)
segments(x0 = data5j1$time, x1 = data5j1$time,
         y0 = data5j1$lower_b, y1 = data5j1$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)

plot(x = data5j2$time,
     y = data5j2$estimate, 
     main = "",
     ylim = c(-25, 15), xlim = c(-0.5, 4.5), 
     xaxp = c(0, 4, 4), yaxp = c(-25, 15, 4),
     pch = 16, cex.axis=1.25)
segments(x0 = data5j2$time, x1 = data5j2$time,
         y0 = data5j2$lower_b, y1 = data5j2$upper_b,
         lwd = 1, col = "black") 
abline(h = 0, lty = 2)


mtext("Vote for the Left", line = -0.1, at = 0.25, outer = TRUE, cex = 1.1)
mtext("Rebel Recruitment", line = -0.1, at = 0.75, outer = TRUE, cex = 1.1)

mtext("Low", line = -1.8, at = 0.13, outer = TRUE, cex = 1.1)
mtext("High", line = -1.8, at = 0.38, outer = TRUE, cex = 1.1)
mtext("No", line = -1.8, at = 0.64, outer = TRUE, cex = 1.1)
mtext("Yes", line = -1.8, at = 0.89, outer = TRUE, cex = 1.1)


mtext("Protest", side = 2, line = .3, at = 0.75, 
      outer = TRUE, cex = 1.1)
mtext("Sanctuary", side = 2, line = .3, at = 0.25, 
      outer = TRUE, cex = 1.1)
mtext(1, text = "Years Relative to the Administration of Treatment", line = 1.2,
      at = 0.5, outer = TRUE, cex = 1.1)

dev.off()

################################################################
# Table 2
################################################################

data <- data_final %>%
  select(dv4b.2.1.farc_ckill_up1_r, iv2.p_strat, iv3.s_strat, iv4b.v_strat, 
         cv1.lpop, cv2.rur_p, cv3c.c_lsh1, cv4.nbi, cv5c.1.bd_o, cv6b.1.ce_o_r, cv7b.3.pgf_cd_o_r, 
         cv8.2.no_p_str, cv8.3.no_s_str, cv8.4b.no_v_str, reb_ru)

stargazer(data, type = 'html', out="Table2.htm")

################################################################
# Table 3
################################################################

data4a <- data4a %>%
  mutate(method = "Protest")

data5c <- data5c %>%
  mutate(method = "Sanctuary")

data1b <- data1b %>%
  mutate(method = "VSP")

data4d1 <- data4d1 %>%
  mutate(method = "Protest_Z4_0")

data5e1 <- data5e1 %>%
  mutate(method = "Sanctuary_Z4_0")

data1f1 <- data1f1 %>%
  mutate(method = "VSP_Z4_0")

data4d2 <- data4d2 %>%
  mutate(method = "Protest_Z4_1")

data5e2 <- data5e2 %>%
  mutate(method = "Sanctuary_Z4_1")

data1f2 <- data1f2 %>%
  mutate(method = "VSP_Z4_1")

data4g1 <- data4g1 %>%
  mutate(method = "Protest_Left_Low")

data5h1 <- data5h1 %>%
  mutate(method = "Sanctuary_Left_Low")

data4g2 <- data4g2 %>%
  mutate(method = "Protest_Left_High")

data5h2 <- data5h2 %>%
  mutate(method = "Sanctuary_Left_High")

data4i1 <- data4i1 %>%
  mutate(method = "Protest_NoRecruitment")

data5j1 <- data5j1 %>%
  mutate(method = "Sanctuary_NoRecruitment")

data4i2 <- data4i2 %>%
  mutate(method = "Protest_Recruitment")

data5j2 <- data5j2 %>%
  mutate(method = "Sanctuary_Recruitment")

table <- bind_rows(data4a, data5c, data1b, 
                   data4d1, data5e1, data1f1, data4d2, data5e2, data1f2, 
                   data4g1, data5h1, data4g2, data5h2, 
                   data4i1, data5j1, data4i2, data5j2)

table <- table %>%
  filter(time == 1)

write_xlsx(table, "Table3.xlsx")

################################################################
# End of this part
################################################################
